/*************************************************
 * Author:      yuanhp
 * Email:       yuan_hp@qq.com
 * 实现fft运算，并将数据写到文件
 * 初始数据从 m.dat 中读取 
 * m.dat 数据格式为每行一个数据
 ************************************************/
#include "iostream"
#include "vector"
#include "fstream"

#include <complex> 
#include <cstring> 
#include <cstdio> 
#include <cstdlib> 
#include <cmath> 

using namespace std;

#define PI 3.14159265358979323846
//#define FS 48000   //48000Hz    //采样率
//-------------- 函数定义 ------------------
void FFT( vector < complex<double> > &TD,  vector < complex<double> > &FD);

int main(){
    vector < complex<double> > ad,ft;
    //读取AD数据
    fstream fin("m.dat") ;
    while(!fin.eof()){
        complex<double> t1;
        fin>>t1;
        ad.push_back(t1);
    }
    fin.close();
    FFT(ad,ft);  //进行FFT变换
    //写FFT后的幅度数据到文件
    ofstream fout("fft_amp.dat");  
    int j=0;
    int r=log2(ad.size());
    long n = 1<<r;
    if(fout.is_open()) {
        for(auto i : ft){ 
#ifdef FS
            fout<<(j)*FS/n<<" "<<abs(i)<<endl;  //只把幅度写入了文件
#else 
             fout<<j<<" "<<abs(i)*2/n<<endl;  //只把幅度写入了文件,FFT后幅度变为了原来幅度的N/2倍，所以除以N/2还原幅度
#endif 
            j++;
        }
        cout<<"已完成"<<j<<"点FFT"<<endl;
        fout.close(); 
    }
    return 0;
}
/******************************************************************************
 * 离散FFT变换
 * @param1：等待变化的离散采样点
 * @param2：变换后的结果
 *****************************************************************************/
void FFT( vector < complex<double> > &TD,  vector < complex<double> > &FD)
{
	long	count;				// 付立叶变换点数	
	int		i,j,k;				// 循环变量
	int		bfsize,p;	
	double	angle;				// 角度	
    
    int r = log2(TD.size());
	count = 1 << r;				// 计算付立叶变换点数
	// 分配运算所需存储器
    vector <complex<double>> W(count / 2,0);
    vector <complex<double>> X1(count ,0);
    vector <complex<double>> X2(count ,0);
    vector <complex<double>> X(count ,0);

 
	// 计算加权系数
	for(i = 0; i < count / 2; i++)
	{
		angle = -i * PI * 2 / count;
		W[i] = complex<double> (cos(angle), sin(angle));
	}
 
	// 将时域点写入X1
	//memcpy(X1, TD, sizeof(complex<double>) * count);
    for(int idx=0;idx<count ; ++idx) {
        X1[idx] = TD[idx];
    }
 
// 采用蝶形算法进行快速付立叶变换
	for(k = 0; k < r; k++)//k为蝶形运算的级数
	{
		for(j = 0; j < 1 << k; j++)
		{
			bfsize = 1 << (r-k);//做蝶形运算两点间距离
			for(i = 0; i < bfsize / 2; i++) 
			{
				p = j * bfsize;
				X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
				X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
			}
		}
		X  = X1;
		X1 = X2;
		X2 = X;
	}
	// 重新排序
	for(j = 0; j < count; j++)
	{
		p = 0;
		for(i = 0; i < r; i++)
		{
			if (j&(1<<i))
			{
				p+=1<<(r-i-1);
			}
		}
        FD.push_back(X1[p]);
	}

}



